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Abstract 

In this paper, the challenging problem of joint channel estimation and data detection for multiple-input 
multiple-output orthogonal frequency division multiplexing systems operating in time-frequency dispersive channels 
under unknown background noise is investigated. Based on two different but equivalent signal models, two 
expectation-maximization algorithm-based iterative schemes for joint data detection and channel and noise variance 
estimation are proposed. The first scheme jointly detects data and estimates the channel and noise variance, but the 
computational complexity is high, owing to the simultaneous detection and estimation for all antennas. To reduce 
the computational complexity, a complexity-reduced scheme that is detecting data and estimating channel for only 
one antenna during each iteration and holding the unknown quantities of other antennas to their last values is 
proposed, whose performance only slightly degrades compared to the first scheme. Moreover, both schemes are 
derived as closed-form expressions, and therefore, our proposed schemes are free of exhaustive search. Simulation 
results demonstrate quick convergence of the proposed algorithm, and after convergence, the performance of the 
proposed algorithm is close to that of the optimal channel estimation and data detection case, which assumes full 
training and perfect channel state information. 

Keywords: Multiple-input multiple-output (MIMO); Orthogonal frequency division multiplexing (OFDM); 
Time-frequency (TF) dispersive channels; Unknown noise variance; Expectation-maximization (EM) 



1 Introduction 

Multiple-input multiple-output (MIMO) communication 
[1] can significantly increase the throughput without 
increasing the transmit power and additional bandwidth. 
Orthogonal frequency division multiplexing (OFDM) 
[2] can provide high data rate transmission capabil- 
ity and is robust against multipath (time-dispersive) 
fading channels. MIMO combined with OFDM (MIMO- 
OFDM) [3] has been adopted in various interna- 
tional standards such as 3GPP-LTE, WiMAX, and 
IMT-Advanced. 
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Meanwhile, vehicles with increased speeds, such as 
high-speed cars, subways, and trains which exceed 350 
km/h, play an increasingly important role in peoples' 
lives. 

Consequently, mobility support is widely regarded as 
one of the key features in current and future wireless 
communication systems. High mobility causes the trans- 
mission channel to change rapidly in time, which results in 
frequency dispersion of the channel. For coherent detec- 
tion in MIMO-OFDM systems, channel state information 
(CSI) is indispensable [3]. 

CSI acquisition is particularly challenging in time- 
frequency (TF) dispersive channels because channel 
responses vary sample by sample, and therefore, the 
number of unknown channel parameters in an OFDM 
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symbol period increases significantly (much greater than 
in frequency-nondispersive channels). Furthermore, in 
practical communication scenarios, the knowledge of 
the power of background noise is required to per- 
form many signal processing algorithms, such as chan- 
nel estimation [4] and decoding [5] in MIMO-OFDM 
systems. 

In this paper, joint data detection and channel and 
noise variance estimation for MIMO-OFDM systems 
operating in TF dispersive channels under unknown 
background noise are investigated. We employ the 
expectation-maximization (EM) algorithm [6,7], which 
is an iterative numerical method employed to com- 
pute the maximum likelihood (ML) estimates, to 
develop an iterative algorithm to solve this challenging 
problem. 

For MIMO systems, the literature along these lines can 
be categorized as follows: 

EM for channel estimation and data detection 
assuming the noise variance is known: EM-based joint 
channel estimation and data detection algorithms 
in time-nondispersive and frequency-nondispersive 
channels (TnDFnD channels) are proposed in [8-10], 
and in time-dispersive and frequency-nondispersive 
channels (TDFnD channels) are proposed in [11-13], 
respectively. However, the maximization step (M-step) 
for data detection proposed in these papers is not 
obtained as a closed-form solution, and therefore, a 
brute-force searching over all of the possibilities is 
required. 

EM for channel and noise variance estimation: In 
TnDFnD channels, EM-based joint channel and noise 
variance estimation algorithms are proposed in [14-16]. 
However, data detection is obtained by an extra ML esti- 
mator and a maximizing a posteriori probabilities (APP) 
detector in [14,15], respectively. 

In [16], a full training sequence is adopted to 
perform the proposed EM algorithm, and there- 
fore, no data detection is addressed. In TDFnD 
channels, EM-based joint channel and noise vari- 
ance estimation algorithms are proposed in [17-19]. 
However, data detection is not addressed in these 
papers. 

EM for data detection and noise variance estimation: 
In TnDFnD channels, an EM-based joint data detection 
and noise variance estimation algorithm is proposed in 
[20]. However, the channel estimate is only obtained by 
pilot symbols and is not included in the EM updating 
process. 

EM only for data detection assuming the noise vari- 
ance is known: In TnDFnD channels, EM-based data 
detection algorithms are proposed in [21,22]. How- 
ever, channel estimation is not addressed in [21], and 
the channel knowledge is assumed ideally known at 



the receiver in [22]. In TF channels, an EM-based 
data detection algorithm is proposed in [23] to solve 
a maximum a posteriori probability (MAP) detection 
problem. However, the data estimate is not given by 
a closed form, and therefore, the exhaustive search is 
required. 

EM only for channel estimation assuming the noise vari- 
ance is known: In TnDFnD channels, EM-based channel 
estimation algorithms are proposed in [24-28]. How- 
ever, the data estimates are obtained by extra MAP 
estimators in [24-26] and APP estimators in [27,28], 
respectively. In TDFnD channels, EM-based channel 
estimation algorithms are proposed in [29-32]. How- 
ever, the data estimates are obtained by an extra 
BI-GDFE detector in [29], a minimum mean-squared 
error (MMSE) detector in [30], a trellises approach in 
[31], respectively, and data detection is not addressed 
in [32]. 

In this paper, based on two different but equiva- 
lent signal models, two EM algorithm-based iterative 
schemes which integrate data detection and channel 
and noise variance estimation are proposed in a con- 
sistent way so as to iteratively improve the system 
performance. 

The first scheme jointly detects data and estimates 
the channel and noise variance, but the computational 
complexity is high, owing to the simultaneous detec- 
tion and estimation for all antennas. To reduce the 
computational complexity of the first scheme, another 
scheme that performs data detection and channel esti- 
mation for only one antenna during each iteration and 
holding the unknown quantities of other antennas to their 
last values is proposed, whose performance only slightly 
degrades compared to the first scheme. Furthermore, 
the estimates of data, channel, and noise variance are 
all obtained as closed-form results, and therefore, the 
proposed schemes are free of exhaustive search. Sim- 
ulation results demonstrate quick convergence of the 
proposed algorithm, and after convergence, the per- 
formance of the proposed iterative algorithm is close 
to that of the optimal channel estimation and data 
detection case, which assumes full training and perfect 
CSI. 

The remainder of this paper is organized as follows. The 
system model for MIMO-OFDM systems operating in TF 
dispersive channels under unknown background noise is 
introduced in Section 2. 

In Section 3, an EM-based scheme for joint data 
detection and channel and noise variance estimation 
is proposed. In Section 4, a reduced complexity EM- 
based scheme is proposed. Section 5 gives some simu- 
lation results that demonstrate the effectiveness of the 
proposed schemes. Finally, conclusions are drawn in 
Section 6. 
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Notation: Matrices and vectors are represented by 
boldface uppercase and lowercase letters, respectively. 

A hat over a variable (e.g., x) indicates an estimate 
of the variable. E{-} denotes the expectation. Super- 
scripts [-] T , and [-] H denote the transpose, the 
matrix inversion, and the Hermitian operations, respec- 
tively. Ijv is an identity matrix with dimension N. diag{x} 
and Blkdiag{-} stand for the diagonal matrix with vec- 
tor x on its diagonal and the block diagonal concate- 
nation of input arguments, respectively. The symbol © 
denotes convolution, and <g> stands for the Kronecker 
product. Tr{X} and |X| are the trace and the determi- 
nant of a square matrix X, respectively. ?fi{-} is the real 
part of the element in the bracket. < • >k denotes 
the mode K operation. The matrix F is the normal- 
ized fast Fourier transform (FFT) matrix with [F] m> „ = 

1 n —i1nmnlN 



xj, = [x p (0),x p (l),...,x p (N p - l)] r denote the data and 
pilot vectors, respectively. 

A cyclic prefix (CP) with length AT cp larger than that 
of the longest channel response is inserted at the begin- 
ning of each OFDM symbol to prevent intersymbol 
interference. 

2.2 TF dispersive channels under unknown background 
noise model 

At the receive antenna /', assuming perfect timing and fre- 
quency synchronization are achieved, the nth sample of 
the received signal is given by: 



N T -1 

y'(«) = J2 #''(«> 0 © 5 '(«) + "A")- ( 2 ) 



2 System model 

2.1 Transmitted MIMO-OFDM systems with scattered 
pilots 

We consider a MIMO-OFDM system with Nt transmit 
and Nr receive antennas. For the z'th transmit antenna, the 
time domain signals' = [s'(0),s'(l), ...,s l (N — l)] r is gen- 
erated by taking the AT-point inverse FFT of the source 
signal in the frequency domain x' = [x'(0),x'(l), ...,x'(N — 
l)] r ass i = E H x i . 
In general, the elements of x' can be categorized into: 



where h )l («, /) is the TF dispersive channel of the Ith path 
with length L at time n, associated with the ith transmit 
antenna and the /th receive antenna, and w'in) denotes the 
unknown background noise and is assumed to obey com- 
plex Gaussian distribution with zero mean and unknown 
variance a 2 , which is assumed to be the same across all 
receive antennas. 

After discarding the CP and stacking all N samples, the 
received signal for a whole OFDM symbol at the receive 
antenna /' can be expressed in a vector form as: 



x'(m) 



x' d (m) V m €l l d 
x l p (m) V m € P p 



(1) N T - 1 

y> = J2 H ; V+w^', 

(=0 



(3) 



where I l d is the index set of subcarriers allocated for data 
symbols (with elements), and 7^, is the index set of sub- 
carriers allocated for pilot symbols (with N p elements), 
respectively. Notice that N = + N p . From (1), we have 
x' = E^x^ + E p x' p , where E l d and E' denote the matri- 
ces collecting columns of Ijv corresponding to I' d and I l p , 
respectively, and x' d = [x' d (0),x' d (l), ...,x' d (N d - l)] r and 



where y> = \y> (0), y> (1), ...,y>(N - l)] r and = 
[w>(0),w>(l),...,w>(N - l)] r denote the received signal 
at the receive antenna /' and the corresponding noise, 
respectively. 

\V l represents the corresponding TF dispersive channel 
matrix and is expressed as: 



#'(0,0) 



#•'(0,1-1). 



H' 1 = 



y\L -1,1-1) W{L - 1,L - 2) . . . W\L - 1, 0) 



#'(0,1) 



0... 



h' l (N -1,1-1) h' l (N - 1,L - 2) . . . y l {N - 1,0) _ 



(4) 
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It is observed from (4) that the number of unknowns 
in W is NL, which is much larger than the num- 
ber of received samples. Therefore, direct estimation of 
H>' is almost impossible (i.e., this will give rise to the 
identifiability problem). 

To overcome this problem, in this paper, a parsimo- 
nious (low-dimensional) representation of h! 1 («, I) using 
the basis expansion model (BEM) [33,34] is adopted, i.e., 
using an expansion with respect to time n of each path I of 
h! l (n, I) into a basis {b„ >q }® =0 as: 

Q 

h ii (n,l) = J2C b r>.<i> ( 5 ) 

q=0 

where P^ is the qth BEM coefficient of the ith path asso- 
ciated with the channel between the ith transmit antenna 
and the ;'th receive antenna; b n>q is the basis that cap- 
tures channel time variations, and Q + 1 is the number of 
the basis. BEM is motivated by the observation that the 
temporal («) variation of h(n, i) is usually rather smooth 
due to the channel's limited Doppler spread and therefore 
{b n>q }q =0 can be chosen as a small set (i.e., Q <?C N) of 
smooth functions. 

Below, two equivalent expressions for the received sig- 
nal will be derived, from which closed-form solution for 
data detection and channel estimation can be obtained, as 
will be shown in the following sections. 

Notice that (3) can be rewritten as: 

N T -1 

y'= G[sW 7 W, (6) 

(=0 

where G[s'] = [diag{s^ 0 }, diag{s* s ^ diagfs^^}] with 
s l csl representing cyclically shifts (cs) s' by i posi- 
tions and h" = [(h£) r ,(hf) r ,...,(hf_ 1 ) r ] r with hf = 
[h' i (0,l),h' i (l,l),...,h' i (N - l,l)] T . (6) can be put into a 
more compact form as: 

y = G[s]h'+V, (7) 

where G[s] = [G[s°] , Gts 1 ] , G[s w i~i]] andb/ = [(h/°) T , 
(b/ 1 ) r , (hi^T-Vff . Using ( 5 ) j h f can b e expressed in 
a vector form as 

hf = Bjf?f, (8) 

where B = [b 0 , bi, bq] withb ? = [b 0iq , bi, q , b N -i q ] T 

and P'l = Wqj, Pi r -, P'qj] 7 ' ■ Substituting (8) into (7), we 
obtain: 

y> = G[s](1 NtL ®B)F+w>, (9) 



where 

fV = [ (j 8/0)T > (j8 ;l ) r > _ ( p(N T -l) )T] T with pi i = [( ^ )T> 

(Pi) T , (Pi_i) T ] T - By stacking the received signals from 
all Nr receive antennas into a single vector using (9) 
and (3), two equivalent expressions of the received sig- 
nal which explicitly show the dependence of the unknown 
BEM coefficient and unknown signal can be obtained, 
respectively, as: 

y=S[s]j8+w (10a) 
= e[j8]s + w, (10b) 

where S [s] = 1 Nr <g> (G[s] (1 NtL ® B)), p = [(P°f, (PV, 
...,{P Nr - 1 ) t ] t , e[ J 8] = [H°,H 1 ,...,H^- 1 ] withH'' = 
[(H oyr (Hi/jH . . . , (H (N R -i )i) H ] H > y = [(y0) r (y i } r _ _ 

iy N R -i ) T ] T i s = [(s o } r (s i ) r | _ _ , iis N T -i )T] t and w = 

[(w°) r , (w 1 ) T , (w Mr - 1 ) t ] t . Notice that S [s] repre- 
sents a function of s and can be reconstructed by s through 
(6), (7), (8) and (9). Similarly, ©[/?] represents a func- 
tion of P and can be reconstructed by P through (3), (4) 
and (5). 

3 Iterative data detection and channel and noise 
variance estimation 

The ML solution of all unknown quantities in (10), i.e., 
s, P, and cr 2 of w, involves multidimensional searches 
that pose prohibitively high computational complexity. In 
this and the next sections, the EM algorithm is employed 
to iteratively compute the ML estimates, with the differ- 
ent accuracy versus complexity trade-offs, respectively. As 
will be seen, our proposed schemes provide not only com- 
putationally affordable but also closed-form solutions that 
are free of exhaustive search. 

Using the EM terminology, we take y as the incomplete 
data, P as the unobservable or missing data, and (cr 2 , s) 
as parameters of interest. The iterative algorithm includes 
two steps (the E-step and the M-step) at each iteration. 
In the E-step, an expectation is taken with respect to P 
conditional on the observed data y and the previous esti- 
mates of (cr 2 , s), and an objective function depending only 
on (cr 2 , s) is obtained. In the M-step, through maximizing 
the function obtained in the E-step, the effect of channel 
can be compensated, and the current updated estimates of 
(cr 2 , s) can be obtained. 
The two steps at the Ath iteration are detailed as follows: 
E-step: compute Q(cr 2 ,s|ff 2 _ 1 ,s j(: _])=E{log/(yj8|cr 2 ,s)|y, 

M-step: solve (a£, s k ) = argmax^ Q(er 2 , s\cr^_ v s^-i). 

Note that conditioned upon y, the only unknown or 
random component in the complete data (y, P) is P; the 
expectation is taken with respect to the conditional proba- 
bility density function f(P |y, a^_ v s^-i), while (cr^ , s^) are 
the estimates of a 2 and s at the kth iteration. 
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More specifically, for the E-step: using Bayes's rule, we 
have: 



f(y,p\a 2 ,s)=f( Y \p,cr 2 ,s)f(P), 



(11) 



is shown in Appendix 1 that the conditional mean and 
covariance matrix are approximately given by: 



(17) 



where the fact that /? is independent of s and a 2 has been 
used. From (11), the function Q(ct 2 , s|cr y 2 _ 1 , s^-i) in the 
E-step can be expressed as: 

Q(<r 2 ,8|o£_i,S*-i) = E{log/(y| J 8,a 2 ,s)|y,cr 2 _ 1 ,s,t_ 1 } 
+ E{log/(/J)|y,cr 2 _ 1 , h-ih 

(12) 

where the second term can be ignored in the following 
derivations, since it is not a function of parameters of 
interest, i.e., not a function of (a 2 , s) and therefore will 
not affect the following M-step. Using (10a), the likelihood 
function /(y|j8, a 2 , s) is obtained as: 



/(yl/W,s) 



n=^ 1 (E if [s M ]S[s w ]) . 



(18) 



■ < n(X 2)N R N 

x exp(--l(y-E[s]/J) H (y- 



H[s]0)J. 

(13) 



Substituting (13) into (12), we have: 

Q{a 2 ,s\a 2 _ v s k -i) oc -NN R \og(na 2 ) 

~ ^2 (y H y-2K{y H S[s]E{^|y,cr 2 _ 1 ,s it _ 1 }} (14) 
+E{/J H E H [s]E[s]/?|y,cr 2 _ p s /< _ 1 }). 

Notice that the following equation holds true for any 
matrix A and vector a with compatible dimension: 



a h A h Aa = Tr{A H Aaa H }. 



(15) 



Define the conditional mean of /? in (14) as ft k = 
E{/?|y, v k _ v s k -i}, and using (15), we obtain: 

Q(<r 2 ,«|o£_i,S*-i) oc -ATAT^logCTrff 2 ) 

--^(^Y-m^ZWh] (16) 
+Tr{E H [s]E[s] (T k + hp")})> 



M-step: in this step, we aim to maximize 
Q(ff 2 ,s|(T i 2 _ 1 ,s J t_i) with respect to a 2 and s. Differenti- 
ating (16) with respect to s and setting the result to zero, 
neglecting those irrelevant terms we have: 

dQ(a 2 , S \a 2 _ v s k -i) 
ds 

x a A j2i1t{y H S[s]^)-Tr{S // [s] S[s] (T*+^")}). 

(19) 

It is noted that since (19) depends on s in an 
implicit way, direct maximization of (19) with respect 
to s is difficult since multidimensional search is 
required. In what follows, an alternative expression for 
Q(cr 2 , s\a k _ v s k -i) will be derived from which a closed- 
form solution for the maximizing value of s can be 
obtained. Since Y* is a N T N R (Q + 1)1 x N T N R (Q + 
1)L Hermitian matrix, based on eigen-decomposition, 
we have f k = ^* (Q+1)I " 1 ^k^m.k^m.k' where 
k m>k is the rath eigenvalue of T k , and fi mk is the 
rath eigenvector, associated with X m , k . Substituting the 
eigendecomposition on Y^ into (19) and using the 
two equivalent equations derived in (10a) and (10b), 
we have: 

29t{y"E [s]p k ] - Tr{E"[s] E[s] (f * + faff)} 

N T N R (Q+1)L-1 

= y"&lP k ]H + s H & H lP k ] y-J2 ^m,kS H ® H lli m ,k] 



m=0 



x ©[^ mtk ]s-s H e H [p k ]©[p k ]s. 



(20) 



Notice that E [•] and ©[•] defined in (10a) and (10b) are 
not only applicable to s and /? but also applicable to any 
vectors with compatible dimension. 

Since (20) is a quadratic form of s, by setting the first 
derivative of (20) with respect to s to zero, the kth signal 
estimate is then given by: 



(N T N R (Q+1)L-1 
m=0 



H=\ X m,k® H [^m,k\®^m,k\ 



where f k = E{G8-/? /< )Q3-/J jt ) H |y, a 2 _ v s^-i} represents 
the corresponding conditional covariance matrix of /?. It 



+& H [hmh] ® H [hh- (2D 
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Note that = [(s° k ) T , (s^) r ,..., (s^ r 1 ) r ] r . After OFDM Recalling (6) and (3), two alternating but equivalent 
demodulation, the symbol from the z'th transmit antenna expressions for y can be derived as: 



can be obtained as: 



N T -l 



(22) 



Since x l k is discrete, belonging to a symbol constella- 
tion point, it must be quantized to its nearest constellation 
point in each iteration. Consequently, constellation map- 
ping is carried out to obtain the discrete symbol estimate 
as: xjj. = Qantfx^.}, where Qantj-} operation denotes quan- 
tization on the element in the bracket. The data symbol 
estimate is thus obtained by collecting the elements of x^. 
corresponding to I' d . 

Finally, putting s* ( = F^xj., i = 0, 1, ...,N T - 1 into (16) 
and setting the first derivative of Q(er 2 , s\& k _ v s^-i) with 
respect to a 2 to zero, the kth estimate of the unknown 
noise variance can be obtained as 



^ = ^ y - m ^ [h]h) (23) 
+Tr{H"[s*]E[s*] (T k + p k p")}y 



In summary, starting from a suitable initial value, the 
proposed iterative EM-based scheme alternates among 
the explicitly closed-form results (17), (18), (21), (22), and 
(23) until convergence, i.e., until no significant changes are 
observed in the updates. 

4 A reduced computational complexity scheme 

The computational complexity of the EM-based itera- 
tive scheme proposed in Section 3 is summarized in 
Table 1. Notice that the computational burden mainly 
comes from the joint detection and estimation simultane- 
ous for all transmit antennas. If in each iteration, detec- 
tion, and estimation can be completed one antenna by one 
antenna, the computational burden will be significantly 
reduced. 



Table 1 Computational complexity of the proposed 
scheme in Section 3 



Computation 


Complexity 


Matrix inversion in (1 7) 


0((N T N R (Q+r)L) 3 ) 


Eigendecomposition on (18) 


0({N T N R (Q+l)L) 3 ) 


Matrix inversion in (21) 


0((N r N) 3 ) 



i=0 



N T -1 

Q[p]s + w= H^Js'+w, 

(=0 



(24a) 
(24b) 



where <t>[s l ] = 1 Nr ® G[s'], H[# c ] = [ (H u <) , (H ) , 
( H (Afe-i)()W]« the BEM coefficients associated with the 
channel from the transmit antenna i to Nr receive 
antennas is represented by fi' rc = [(P°') T , (j8 u ) r , 
yito-Diff with p = {{ f Q )T t ^ i) T t „.,(^_ l) r ] r 

The subscript 'rc' is short for 'reduced complexity' 
to distinguish it from the /? ; defined in Section 3, 
which represents the BEM coefficients associated with 
the channel from Nt transmit antennas to the receive 
antenna /'. 

From (24a) and (24b), it is observed that by applying 
the mathematical framework of EM, an alternative way to 
choose the complete data, defined as ^ in this scheme, 
is by decomposing the observed data y into its signal 
components. The complete data is obtained as: 



O[s 0 ] (l NRL ®B)p 0 rc 
*[s'](W8B)j!! c 



Lfcfs^KlAfcLSB^-iJ |_w 



H[/5? c ]s° 





r w o i 


+ 


w 1 




_ w A?'r-l_ 







r w° -| 




+ 


w 1 


1 




w ^-i 



(25a) 
(25b) 



where f = 0[s ; ] (1 NrL <g> B)j8^ c + w< = H[# c ] s'+w', i = 
0, 1, ...,Nt — 1, and w', i = 0, 1, ...,Nt — 1 are circularly 
symmetric and statistically independent Gaussian vectors 
satisfying w = Ylflo' 1 W. 
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Similar to the E-step in Section 3, for the kth itera- 
tion, we need to compute the conditional expectation of 
the log-likelihood function for the complete data ijr. More 
specifically, for the 

E-step: using (25a), the likelihood function can be 
expressed as: 



f(t\P rc ,s) = 



xexp(-(f-Q) H r- l (i/r-Q)y 



(26) 



where /J rc = [(p° tc ) T , (# c ) r , (jS^ 1 ) 7 ] 7 , "*V = 
Blkdiag{ ? °CT 2 I WNi? , s^Hmn,,, .-, ^-V 2 ^} with 
E{w'(w') H } = g l a 2 lMN R and the g-"s being arbitrary non- 
negative and real-valued scalars satisfying Y^aIq 1 S"' = 1> 
Q = [($[s°] (l NsL <g> B)^° c ) r ,(<I)[s 1 ] (I^i ® B)/?J c ) r , 
($[ s ^r-l] ^ B)^ r " 1 ) r ] r . Notice that in this 

scheme, we take (/3 rc ,s) as parameters of interest. 
Using (26) and neglecting those irrelevant terms, 
E{log/(^|/? rc ,s)|y, Pick-vH-i] can be expressed as: 



E{logf(t\p rc ,s)\y,p ICik _ v s k -i} 



xcxQ H r- l f k + f"r- 1 Q -Q H r- 1 6 , 



(27) 



where ijr k , the conditional mean of can be derived as: 



tk = Mf\Y,Prc,k-l>h-l} 

= Q k _, + r f i5 H (i5r f i5 H r\y- 



(28) 



where 13 = [Inn r > Inn r > — > Inn r ] is a matrix with dimen- 
sion NNr x NNrNt that connects y and i/r as y = 
YlfJo 1 = Substituting the corresponding com- 
ponents into the right-hand side of (28), after some 
manipulations we obtain: 



Substituting (29) into (28), finally we obtain: 



E{log/(^|/? rc , s)|y, P rC:k _ v s k -i) 

N T -1 



(**-fl>[s'] (iAfciOB)^), 



(30) 



(=0 



where 



t k = <t>[s k _ l ](l Ns L®B)P [C ,k-i 

Nt-i (31) 

+ s'(y- E *[»U]( 1[ A6 t i®B)^_ 1 ). 



(=0 



It is noted from (30) that in the following M-step, 
the maximization of E{log/\^|/? rc , s)|y, j8 rc> ^_i,s^_i} with 
respect to /? rc and s is equivalent to the minimization 
of each of the single terms in (30), i.e., minimization of 

(**-4>[s'] (lN R L®B)Pi c ) H (f k -<t>ls i ] dAfciBB^with 
respect to /?J. C and s' for each i, separately. 

Notice that the multidimensional minimization for each 
of the terms in (30) still remains a formidable task. To 
solve this problem, substituting (24a) and (24a) into (31), 
we obtain: 



fk 



where 



S i <!>ls i ](l NsL ®B)Pi c + X i k 



(32a) 
(32b) 



4=^*k-l\( l N S L®B)h c ,k-l 
I N T -1 

N T -1 \ 
- E ^H-l^ R L®m'rc,k-X 

(33) 



i=0 



O^KlAki ® B)^ a _! + c; 0 (y - ES) 1 f^-Jd^i ® B)^-l) 
Ots^Jd^i ® B)^_! + c; 1 (y - ES)" 1 <*> [s[_J(W <8 B)^ a -i) 



(29) 
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Set g l = 1, for the ith transmit antenna we have: 

f\ = 4>[s'](lAfc£ <8> B)# c + xi (34a) 
= H[#Js' + x£. (34b) 

Recalling that S - ' = an d by usm g (31). for the 

transmit antennas {g}^^,- we have: 



Using (34) and (35), (30) can be decomposed into Nt 
terms, each of which can be solved as follows: 
M-step: for the ith transmit antenna, we have: 

[Kck>K] = ar gmin j^-^Mto,! ® B)# c )" 

x (^-^[s'laAwOB)^)), 

(36) 



and for the transmit antennas Lg^Jn^f' we have: 
[^Kyt'^f] = [/*rc,it-l> sf_i] • 



(37) 



Therefore, the proposed iterative scheme starts from 
k = 0, 1, 2, ... and during the /cth iteration, i is set as i = 
< k >n t - It can be seen that we have split the estimation 
and detection problem for the MIMO case of Section 3 
into estimation and detection problem for Nt single-input 
and multiple-output (SIMO) cases, where, during each 
iteration, parameters and data from only one transmit 
antenna are estimated and detected. Note that x.\ given 
in (33) is a disturbance term that accounts for the back- 
ground noise and residual interference after the Ath iter- 
ation, where the interference is linearly related to the sig- 
nals of all transmit antennas. Then, assuming the interfer- 
ence is i.i.d with zero mean, from the central limit theorem 
[35], it can be seen that the entries of x' k are nearly Gaus- 
sian distributed with zero mean and some variance er 2 , , . 
Under the above assumption, it turns out that the mini- 
mization problem in (36) is equivalent to the ML estima- 
tion of f}' rc , s l and the unknown variance cr 2 ,- starting from 

~ i 

the observation f . Comparing (34a) and (34b) with (10a) 
and (10b), it is easy to see that the same EM procedure 
proposed in Section 3 can be directly adopted to solve 
the optimization problem of (36), with details shown in 
Appendix 2. 

The computationally feasible EM scheme is 
summarized as follows: 



Iteration until convergence: 
for the /cth iteration do 

Set the transmit antenna as i = < k >m t 

1. Updating ijr 1 : 

N T -1 

g=o,g& 

1. Updating ji l rc , t\ c , s l and a 2 t until convergence 



(35) (assume the convergence occurs at the £th iteration): 



for the /th iteration do 
Updating j} l rc and Y l rc : 

Kk-U = * x Vu-i ((B^G^s^.J 

where s[_ 1>0 = s\_ v a 2 .^ = oj^, +\ = 

KitY, (*JV, .... (fT~ U ) T V and B z = i£ <g> B. 

Then, P^-ij can be obtained as 0 rCj £-i,/ = 

[()8ra-i,/) r .()8 r a-u) T ,-.(i8 r a- W ) r ] r and 

"^rc/<-i,/ can be obtained as Tf rc ;t_i/ = 
Blkdiagff t^_ i; , Y^_y}. 

Updating s l : 

(N R (Q+1)L-1 
E Kcm.k-U^lKcm.k-u] 
m=0 

X H[H^_ W ]+H"^_ U ] 

wherp f' _ v JV«(Q+l)i-l , / « 

wnere I TcJc _ U - Lm=0 \c,m,k-ljl l rc,m,k-U 

(i L [ cm j ( _i [) H represents the eigendecomposition of 

1 rc,/<-l,/- 



Updating a 2 , : 



V 1 



((^)"^-2((^)"0[s'_ u 



x (\ NrL ® B)p [Cik _ u \ + Tr {(Iam® B) 
x 3>"[s[_ 1(/ ] 0[s[_ 1(/ ] (I^OB) 

X (^rc/t-U + /W-W^W-l,/)^)}) • 



H 



end for 
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3. Set [p rCik ,s l k ,a^] = W r ck-l,£'S k -l,£' a x',k-l,C ] f 0r 
the transmit antenna i. 

Set [f rcM , sf, b\ t £ = itftf-L 4-v f° r the 

transmit antennas {jj}^o^«- 
end for 



The computational complexity of the proposed iter- 
ative EM-based scheme with reduced complexity is 
summarized in Table 2. 

Note that compared to Table 1, the computational com- 
plexity of the proposed iterative EM-based scheme with 
reduced complexity is significantly lower than that of 
the EM-based scheme proposed in Section 3. However, 
this significant computational complexity reduction is not 
obtained without price. As will be shown in Section 5, 
there is a minor performance degradation compared to 
the EM-based scheme proposed in Section 3. This per- 
formance degradation is due mainly to two reasons. First, 
the disturbance term in (34a) and (34b) contains the back- 
ground noise as well as the residual interference from 
other transmit antennas, whereas in (10a) and (10b), only 
the background noise is contained. Second, the separate 
estimation and detection for each antenna is seen as a 
suboptimal estimation and detection method compared 
to the joint estimation and detection for all antennas, 
which is optimal in the sense of estimation and detection 
theory [36]. 

4.1 Initialization 

The EM algorithm is guaranteed to obtain at least a local 
maximum after convergence [6,7]. 

To provide an initial value, a least square (LS) algorithm 
based on pilot symbols is utilized to provide a good initial 
estimate which will be demonstrated in the simulations. 
Recalling (1), (10a), and (10b), we have: 



y = S [St p x p ]p + Blfi] Sl d x d + w, 



(38) 



where Q p = Blkdiag{F H E£, V M E p , F^E^" 1 }, Q, d 
Blkdiag{F"E°, F"£i, ^ N /~\ x p = [(x°) r , {^f, 

(xf -yr, and x d = [(x°) r ,(xi) r ,...,(xf-yF. By 



treating the term containing as interference, the LS 
estimate of /? is obtained as: 

P 0 = {^"[Qp^ElQpXp])' 1 E H [a p x p ]y. (39) 

Substituting (39) into (10b), the initial signal detection 
is obtained as: 



so = (In t ® F H )Qant 



(i A , r ®F)(e H [^ 0 ]e[ J 8 0 ]) 1 



(40) 

Finally, for the initial variances CTq and {o^ ( 0 }jlo _1 > tne y 
are all set to 0. 



5 Simulation results and discussions 

In this section, the performance of the proposed algo- 
rithm is demonstrated by Monte Carlo simulations. In 
the simulations, transmit and receive antennas are set as 
Nt = Nr = 2, each OFDM symbol has 64 subcarriers 
(AT = 64) and communicates over a bandwidth of 20 MHz. 
The sampling interval T s is thus 50 ns. The length of the 
CP is AT cp = 8. 

The normalized maximal Doppler shift is set as Nf d T s = 
0.075 and 0.15, respectively, where f d represents the max- 
imum Doppler frequency. 

The channel has three taps (L = 3) with an exponen- 
tial power delay profile, namely erf = exp(— kI)((1 — 
exp(-K))/(l-exp(-KL))),/ = 0, 1, L - 1 with k = 1/3. 
In typical communication scenarios, only a few signifi- 
cant paths dominate the effect of the wireless channel 
[4]. Therefore, L = 3 is a reasonable setting. Each 
tap coefficient follows a complex Gaussian distribution. 
The data are modulated by quadrature phase shift key- 
ing (QPSK) and 16 quadrature amplitude modulation 
(16 QAM), respectively, with unit power. The pilot clus- 
ter follows the structure in [37], and more specifically, 
seven pilot clusters are used for each transmit antenna. 
The clusters are equal-spaced among subcarriers, and in 
each cluster, one nonzero pilot is guarded by one zero pilot 
on each side. The nonzero pilots are generated as zero- 
mean complex Gaussian random variables with power 
three times that of data symbols. Furthermore, the gen- 
eralized complex exponential BEM (GCE-BEM) [34] is 
adopted. 



Table 2 Computational complexity of the proposed 
scheme in Section 4 



Computation 


Complexity 


Matrix inversion in (56) 


O(((Q+1)0 3 ) 


Eigendecomposition on (57) 


O(((Q+1)0 3 ) 


Matrix inversion in (53) 


0(W 3 ) 



5.1 Convergence of the proposed schemes 

Figures 1, 2, and 3 present the convergence perfor- 
mance of the proposed EM-based scheme in Section 3 
(marked as scheme 1) and the proposed EM-based 
scheme in Section 4 (marked as scheme 2) with signal- 
to-noise ratio (SNR) equal to 10, 20, and 30 dB. It 
can be seen that both the mean-square error (MSE) 
and bit error rate (BER) improve significantly in the 
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Scheme 2 (SNR=10dB) 
Scheme 1 (SNR=10dB) 
Full training (SNR=10dB) 
Cramer Rao bound (SNR-lOdB) 
Scheme 2 (SNR=20dB) 
Scheme 1 (SNR-20dB) 
Full training (SNR=2()dB) 
Cramer Rao bound (SNR=20dB) 
Scheme 2 (SNR=30dB) 
Scheme 1 (SNR-30dB) 
Full training (SNR=30dB) 
Cramer Rao bound (SNR=30dB) 



7 8 9 10 
Number of iterations 

Figure 1 Convergence of channel estimation, QPSK, 

NfyTs = 0.075. Convergence performance of the proposed 
EM-based scheme in Section 3 (marked as scheme 1) and the 
proposed EM-based scheme in Section 4 (marked as scheme 2) with 
signal-to-noise ratio (SNR) equal to 1 0, 20, and 30 dB. 




Scheme 2 (SNR=10dB) 
Scheme 1 (SNR=10dB) 
Perfect CSI (SNR=10dB) 
Scheme 2 (SNR=20dB) 
Scheme 1 (SNR=20dB) 
Perfect CSI (SNR=20dB) 
Scheme 2 (SNR=30dB) 
Scheme 1 (SNR=30dB) 
Perfect CSI (SNR=30dB) 



H H 



-+- 



•0 0 © 0 0 6 © 0- 



v — v — v — v — v — v — v — v- 



7 8 9 10 
Number of iterations 



11 12 13 14 15 



Figure 3 Convergence of data detection, QPSK,Nf d T s = 0.15. 

Convergence performance of the proposed EM-based scheme in 
Section 3 (marked as scheme 1) and the proposed EM-based scheme 
in Section 4 (marked as scheme 2) with signal-to-noise ratio (SNR) 
equal to 10, 20, and 30 dB. 



first few iterations and converge to stable values within 
eight iterations. Channel estimation with full training and 
data detection with perfect CSI are shown for compar- 
ison. Furthermore, according to [38], the Cramer-Rao 
bound is also shown for comparison. It can be seen 




Scheme 2 (SNR=10dB) 
Scheme 1 (SNR=10dB) 
Perfect CSI (SNR=10dB) 
Scheme 2 (SNR=20dB) 
Scheme 1 (SNR=20dB) 
Perfect CSI (SNR=20dB) 
Scheme 2 (SNR=30dB) 
Scheme 1 (SNR=30dB) 
Perfect CSI (SNR=30dB) 



-e — e — e — e — e — o — e — e — e- 
— v — v — v — v — v — v — v — v- 



Number of iterations 

Figure 2 Convergence of data detection, QPSK,NfdT s = 0.075. 

Convergence performance of the proposed EM-based scheme in 
Section 3 (marked as scheme 1) and the proposed EM-based scheme 
in Section 4 (marked as scheme 2) with signal-to-noise ratio (SNR) 
equal to 10, 20, and 30 dB. 



from Figure 1 that after convergence, the channel esti- 
mation performance of both schemes greatly improve 
that of the initial estimation (marked as iteration = 0), 
which indicates the ability of the proposed algorithm to 
cancel the interference from unknown data to channel 
estimation through iterations. The channel estimation 
performance of scheme 1 is very close to that of the 
Cramer-Rao bound and the full training case. The chan- 
nel estimation performance of scheme 2 suffers a minor 
performance degradation compared to that of the scheme 
1, which is the price we have to pay for the reduced com- 
putational complexity. Similar results can be observed for 
the performance of data detection in Figures 2 and 3, 
which indicates that the updated channel estimate can 
in turn greatly improve the data detection through iter- 
ations. Similar convergence results are also observed for 
the 16 QAM case, and figures are not presented here due 
to space limitations. 

5.2 Performance of the proposed schemes 

Figures 4, 5, and 6 show the MSE and BER performance 
achieved by the proposed iterative algorithm versus SNRs. 
It can be seen from Figure 4 that the performances of the 
proposed schemes 1 and 2 both perform much better than 
that of the initial value and close to that of the Cramer-Rao 
bound and the full training case after convergence. 

Similarly, it can be seen from Figure 5 that for the case 
where Nf^T s = 0.075, the BER performance of the pro- 
posed iterative algorithm is very close to that of the ideal 
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Initial value (ileralion=0) 
Scheme 2 (ileralion=l) 
Scheme 1 (ileralion=l) 
Scheme 2 (ileralion=8) 
Scheme 1 (ilcralion=8 
Full Iraining 
Cramer Rao bound 



Figure 4 Performance of channel estimation, QPSK, 

NfjTs = 0.075. MSE performance achieved by the proposed 
iterative algorithm versus SNRs. 





Figure 6 Performance of data detection, QPSK and 1 6 QAM, 
NfaT s = 0.15. BER performance achieved by the proposed iterative 
algorithm versus SNRs. 



case which assumes perfect CSI after convergence. For 
the severe case where NfdT s = 0.15, it can be seen from 
Figure 6 that the proposed iterative algorithm can still deal 
with such a highly TF dispersive channel and performs 
well. Moreover, from Figures 5 and 6, it can be seen that 
for signals with both amplitude and phase variations such 
as 16 QAM, the proposed algorithm also performs well. 



Finally, we investigate how the proposed schemes are 
affected by different channel lengths. A severe case 
where the channel length is equivalent to the number 
of embedded pilots (marked as case 2) is shown in 
Figure 7. As can be seen from the figure, compared 
to the originally-presented case where the channel 
length is 3 (marked as case 1), there is an obvious 




16QAM(inilial value, ileraliou=0) 
16QAM(scheme 2. iteration=l) 
16QAM(schcmc 1. itcration^l) 
16QAM(schcmc 2. itcration=8) 
16QAM(schcmc 1. itcration=8) 
16QAM(perfect CSI) 
QPSK(initial value, iteration^)) 
QPSK(scheme 2. ileralion=l) 
QPSK(schcmc 1, iteration=l) 
QPSK(scheme 2. iteration=8) 
QPSK(scheme 1. ileralion=8) 
QPSK(perfecl CSI) 



1() 



15 

SNR.dB 



Figure 5 Performance of data detection, QPSK and 16 QAM, 

NfdT s = 0.075. BER performance achieved by the proposed iterative 
algorithm versus SNRs. 




Case 2 (initial value, iteration=0) 
Case 2 (scheme 2. ileration=8 
Case 2 (scheme 1. iteration=8 
Case 1 (initial value. iteralion=0) 
Case 1 (scheme 2. iteration^) 
Case 1 (scheme I, iteration=8) 
Perfect CSI 



10 



15 
SNR.dB 



Figure 7 Performance of data detection for different channel 
length, QPSK,jV/ d r s = 0.075. BER performance achieved by the 
proposed iterative algorithm versus SNRs. 
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performance degradation of the proposed schemes for 
the severe case 2. The reason can be explained according 
to the estimation theory [36] that when the channel 
length increases, more parameters need to be estimated, 
which leads to a decreased performance. On the con- 
trary, if the channel length decreases, less parameters 
need to be estimated and that leads to an increased 
performance. 

6 Conclusions 

In this paper, two EM-based iterative data detection 
and channel and noise variance estimation schemes 
for MIMO-OFDM systems operating over TF disper- 
sive channels under unknown background noise have 
been proposed. The resulting schemes achieve conver- 
gence in a few iterations and can effectively estimate 
TF dispersive channels and obtain reliable data detec- 
tion under unknown background noise environments. 
The first scheme iteratively detects data and estimates the 
channel and noise variance simultaneously for all anten- 
nas, and moreover, the updating expressions of these esti- 
mates are all derived as closed-form results. Simulation 
results showed that after convergence, the performance of 
the first scheme is very close to that of the optimal case 
which assumes full training and perfect CSI. To reduce 
the computational complexity of the first scheme, another 
EM-based scheme that detecting data and estimating 
channel for only one antenna during each iteration and 
holding the unknown quantities of other antennas to their 
last estimates has been proposed, which is also derived 
as closed-form results. Simulation results showed that its 
performance only slightly degrades compared to the first 
scheme, but the computational complexity is significantly 
reduced. 

Appendices 

Appendix 1 

Derivation of (17) and (18) 

Using Bayes's formula, the conditional pdf of /? is given by: 



/(/8|y»o}t_i»8*_i) 



/( y |/J,cr 2 ,s)/(/?) 



f(y\(r 2 ,s) 



(T 2 =5' jt 2 _ 1 ,S=SA : _l 



(41) 



where the fact that (i is independent of s, and a 2 
has been used. The BEM coefficient /? can be shown to 
be complex Gaussian variable [33] with zero mean and 
covariance matrix R^, that is: 



f(fi) = ttwttjtt-, r expC-jS^RI 1 /?). 



Note that: 

/(y|a 2 ,s) = jf(y\p,cT 2 ,s)f(p)dp. 



(43) 



With f(y\P, a 2 , s) given by (13), putting (13) and (42) into 
(43), we have: 



/(yk 2 ,s) 



\r k \ 



(7T<T 2 ) NN z\Rf)\ 



x expt-^Cy^y-a^ff/^)). (44) 



Substituting (13), (42), and (44) into (41), after some 
manipulations we have: 

i 

/(/J|y»o£_i»8*_i) 



n (Q+W T N R L\Y k \ 



~ -1 



(42) 



x exp(-(/? - f} k )"r k {fi - fi k )), 

(45) 

where 

fi k = (E H [s k - 1 ]E[s k -i] + a 2 _ 1 R/y 1 E H [s k -i]y, 

(46) 

f k = ^(S^fe-i] S + o^R" 1 )- 1 . (47) 

Thus, the pdf f(fi\y,cr k l ,s k -i) is a Gaussian distri- 
bution. In addition, fi k and T k given in (46) and (47), 
respectively, are in fact its conditional mean and covari- 
ance. To show that we have no prior information on /?, we 
take the limit 1 1 R/s 1 1 — »- +oo, which leads to (17) and (18). 
In this paper, we set R^ 1 to zero to show we have no prior 
information for /?. Indeed, there will be a performance 
degradation by assuming R^" 1 to zero. However, this is a 
typical complexity versus performance trade-off. More- 
over, as can be seen from simulation results in Section 
5, even we set R^ 1 to zero, the proposed algorithm also 
performs well, and its performance is acceptable. 

Appendix 2 
Solving (36) 

Comparing (34a) and (34b) with (10a) and (10b), referring 

to Section 3, we take f' k as the incomplete data, f}' IC as the 
unobservable or missing data, and (a 2 it s') as parameters 
of interest. The two steps at the £th iteration are detailed 
as follows: 

E-step: compute Q^.s^^.sj^) = E{log/X^, 
M-step: solve (a 2 ik ,s l k ) = argmax^^Q^.s 4 !^,^^ 
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Note that conditioned upon ijr k , the only unknown or 

random component in the complete data (rjr k , fi\ c ) is f} l IC , 
the expectation is taken with respect to the conditional 

probability density function/ (/? [ c \ ijr k , k _ 1 , s' k _ x ) , while 
(a 2 , ,,s 1 ,) are the estimates of a 2 ,, and s' at the Ath itera- 
tion. More specifically, for the E-step: Using Bayes's rule, 
we obtain: 

/(^ rc |a x V) =/(fil/4.^V)/(# c )- (48) 

Using (48), the function Q(ct 2 ,, s'|cr 2 ^ ,>sl ,) can be 
expressed as: 

Q(<#.«'|o$jt_i.8jLi) 

= E{log/-(^l# c > trj, s')|^{, a 2 ^, s^j} 

+E{iogf^^)if;, c > x 2 , > _ 1 ,s[_ 1 }, 

(49) 

where the second term can be ignored in the following 
derivations, since it is not a function of parameters of 
interest, i.e., not a function of (cr 2 ,, s'). Using (34a), the 

likelihood function f(ijr k \p l rc , er 2 ,, s') is obtained as: 

x' 

x exp(-4"(^ - * [»'] dN R L ® B)/?' c )" 

/'* 

x (**-*[s'] (Iiv K i®B)^ c )). 

(50) 

Substituting (50) into (49) and referring to (14), (15) and 
(16) and Appendix 1, the conditional mean and covariance 
matrix are obtained as: 

Kc,k = ((W <8> B)"*"^] Ots^] (1 NrL ® B))" 1 
xfl^SBy^J^ (51) 

x(1 NrL ®B))-\ (52) 

It is noted that the matrix [s ( ] is of dimension NrN x 
NrNL, and the matrix (1n r l ® B) is of dimension NrNL x 
Nr(Q+1)L; theNn(Q+l)LxNn(Q+l)L matrix inversion 
required in (51) and (52) is only ^ of that needed in (17) 
and (18). 



M-step: using the two equivalent expressions derived in 
(34a) and (34b) and similar to (19), (20), (21) and (22), the 
signal updating equation is obtained as: 

A/«(Q+l)£-l 

~ S 'k = ( A rc, m ,£ HW [/4c,m,*] H [/4c,m,/<] 

+ H"[£ a ]H[£ a ] )~ l n H [^ k ] t\, (53) 

Wnere I rc ^. 2__, m= o A rc,m,kf l rc,m,kWrc,m,k> re P 

resents the eigendecomposition of T rc jt . It is noted that 
compared to (21) where NfN x NtN matrix inversion is 
required, only N x N matrix inversion is needed in (53). 
The symbol detection can thus be obtained after OFDM 
demodulation as 

x[ = Qant{F4}. (54) 

Substituting = ¥ H x k and (50), (51) and (52) into (49) 
and referring to (23), the unknown noise variance for the 
disturbance term x' k can be obtained as: 

+ Tr{(l NsL <g> B) H $ H [s k ]<i>[s k ] (1 NrL ® B) 
x(tc,k + KcAKc,kf)})- 

(55) 

In summary, (51), (52), (53), (54) and (55) solve the 
minimization problem in (36). 

Notice that the computational complexity can be further 
reduced by observing the diagonal structure of both <t> [s l ] 
and (1n r l® B) in (24a). Therefore, (51) and (52) can be fur- 
ther split into Nr sub-matrices, each of which is expressed 
as: 

££t = ((B*)"G ff [i' t _ 1 ]G[a' t _ 1 ]B«)- 1 (B*) H G"[a' t _ 1 ]^, 

(56) 

<k= ^,.- 1 ((B Z ) H G H [sj < _ 1 ]G[sj ( _ 1 ] B^)" 1 , (57) 

where f k = [(f k f, (f k f , (f k f] T . Then, p lck 
is obtained as 0 rcJt = [(/5 ra ) r , (/? ra ) r , .... (j8 rc>(t ) r ] r , 
and T rc j, can be obtained as T rc k = Blkdiag{T^' /( , Y^, 
Y^ -1 ''}. It is noted that the matrix G[s ! ] is of dimen- 
sion N x NL, and the matrix B z = li (g) B is of dimension 
NL x (Q + l)L; the (Q + l)i x (Q + l)i matrix inversion 
required in (56) and (57) is ^ of that needed in (51) 

and (52) and therefore only N * N of that needed in (17) 
and (18). 
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